Introduction
A Partial Differential Equation (from now on, PDE) is an equation where both functions and their partial derivatives may appear.
PDEs are almost always used to study and predict physical phenomena: they usually consists of multiple variables and a close solution almost never exists, so we need to use computer to perform all the calculations.
Assume you want to study some physical phenomena, you cannot feed it directly into the computer: first you have to create a mathematical model of the phenomenon, deriving the euqtions that describe it (modelling phase), then you have to discretize the model so that it can be treated with numerical methods, obtaining an algebraic system (discretization phase) and only then you can feed it into the computer (computation phase).
Each of those three phases introduces an error: let u_{ph} be the real solution (what you can observe in the real world), u_m be the error introduced in the modelling phase, u_n the error introduced in the discretization phase and u_c the error introduced by the finite-precision computation, then
\|u_{ph} - u_c\| \le \underbrace{\|u_{ph} - u_m\|}_\text{modelling error} + \underbrace{\|u_m - u_n\|}_\text{numerical error} + \underbrace{\|u_n - u_c\|}_\text{computation error}
u_c is the solution we effectively get out of all the simulation.
In this document we assume that there is no modelling error.
There exist three types of PDEs: elliptic, parabolic and hyperbolic. While the first type only depends on space, the other two depend also on time. Hyperbolic PDEs will not be considered in this document.
First, an overview about PDEs will be given and then, numerical methods for the approximation of their solution will be analyzed.
Elliptic PDEs
A general elliptic PDE can be expressed as
\begin{cases} Lu = f & \text{in } \Omega \\ \text{Boundary conditions} & \text{in } \partial\Omega \end{cases}
where \Omega \sub \mathbb{R}^d, d \in \mathbb{N}^+ is called the domain of the PDE and \partial\Omega is the boundary of \Omega, f : \Omega \to \mathbb{R} is a given function and u : \Omega \to \mathbb{R} is the unknown function.
Tipically, the boundary condition can be Dirichlet conditions and Neumann conditions:
\begin{cases} u = g & \text{in } \Gamma_D \sube \partial\Omega \\ \frac{\partial u}{\partial n} = q & \text{in } \Gamma_N \sube \partial\Omega \end{cases}
where g : \Gamma_D \to \mathbb{R} is called Dirichlet data and q : \Gamma_N \to \mathbb{R} is called Neumann data. \frac{\partial}{\partial n} denotes the normal derivative (see appendix).
\Gamma_D is called Dirichlet boundary, \Gamma_D is called Neumann boundary. They form a partition of \partial\Omega:
\begin{cases} \Gamma_D \cup \Gamma_N = \partial\Omega \\ \Gamma_D \cap \Gamma_N = \empty \end{cases}
A problem can be a Dirichlet problem (if it has only Dirichlet conditions), Neuman problem (if it has only Neumann conditions) or mixed problem.
L is a differential operator called elliptic operator and it is defined as a sum of three terms (in order, the diffusion term, the advection/convection/transport term and the reaction term):
Lu \overset{\Delta}{=} -\operatorname{div}(\mu \nabla u) + \vec{b} \nabla u + \sigma u
For this reason L is also called advection/diffusion/reaction (ADR) operator.
See appendix for an explanation of the various symbols.
If \mu is constant then -\operatorname{div}(\mu \nabla u) = -\mu \Delta u.
Monodimensional elliptic PDEs
In order to solve an elliptic PDE, the problem needs to be massaged to get to something that can be treated with computers.
flowchart LR
id1[Strong problem]
id2[Weak problem]
id3[Galerkin problem]
id4[Algebraic problem]
id1-->id2-->id3-->id4Getting to the algebraic problem
The strong problem is the problem as it was defined before:
\text{Find $u$ s.t.} \begin{cases} Lu := -(\mu u')' + b u' + \sigma u = f & \text{on } \Omega \\ \text{Boundary conditions} \end{cases}
where all the derivatives are weak derivatives.
The strong problem always has one unique solution.
In order to get to the weak problem we have to multiply both sides of the PDE by v \in V and then integrate them on \Omega:
\begin{align} \int_\Omega -(\mu u')' v + b u' v + \sigma u v &= \int_\Omega f v \\ \int_\Omega \left( (\mu u')' v' + b u' v + \sigma u v \right) - \left[ \mu u' v \right]_{\partial\Gamma} &=\int_\Omega f v \end{align}
where V is defined as follows:
V = \left\{ v : \Omega \to \mathbb R \middle| v \in L^2(\Omega), v' \in L^2(\Omega), v(\Gamma_D) = 0 \right\}
where L^2 is the set of square-integrable functions.
The term in the square brackets has to be evaluated immediately: it is 0 on the Dirichlet boundaries and a real number in the Neumann boundaries. The terms that comes out from the evaluation in the Neumann boundaries (\sum \gamma v(x_D)) must be moved to the right hand side of the equation.
The weak problem is then written as follows:
\text{Find $u \in V$ s.t. } \underbrace{\int_\Omega (\mu u' v' + b u v' + \sigma u v) dx}_{a(u, v)} = \underbrace{\int_Omega f v \, dx+ \sum \gamma v(x_D)}_{F(v)} \qquad \forall v \in V
that is equivalent to
\text{Find $u \in V$ s.t. } a(u, v) = F(v) \qquad \forall v \in V
a is a bilinear form and F is a linear functional.
As the weak problem must be solved \forall v \in V, it has infinitely many solutions. We do not like that, so we restrict the space in which the solution lies to a smaller, finite dimension one.
Let V_h \sub V be a finite-dimension function space, then the problem can be rewritten in the Galerkin form
\text{Find $u_h \in V_h$ s.t. } a(u_h, v_h) = F(v_h) \qquad \forall v_h \in V_h
Since V_h has a finite number of dimensions, then it admits a basis \{\varphi_j\}_{j=1}^{N_h} such that
u_h(x) = \sum_{j=1}^{N_h} u_j \varphi_j(x)
By exploiting the linearity of a and F, then the problem can be rewritten again as multiple simpler equations, one for each element of the basis:
a\left( \sum_{j=1}^{N_h} u_j \varphi_j, \varphi_i \right) = F(\varphi_i) \iff \sum_{j=1}^{N_h} u_j a(\varphi_j, \varphi_i) = F(\varphi_i) \qquad \forall i = 1, 2, \dots, N_h
This can also be rewritten as A \vec{u} = \vec{F} (algebraic form). The algebraic form is a simple (and, usually, megachonk) linear system so you can have fun solving it.
Construction of V_h
Let
\begin{align*} \mathscr{T}_h &= \left\{ k_j : \bigcup_{i=1}^{N_h} k_i = \Omega, k_i \cap k_j = \empty \ \forall i, j \in 1, 2, \dots, N_h, i \ne j \right\} \\ X_h &= \left\{ v_h \in \Omega \to \mathbb{R} \middle| v \in \mathcal{C}^0(\Omega), {v_h}\big|_{k_j} \in \mathbb{P}^g(k_j)\ \forall j \right\} \end{align*}
We split \Omega in N_h h-wide elements.
We want v_h \in V_h to be a piecewise polynomial (i.e. continuous in each element) of degree g on \Omega such that v_h(\Gamma_D) = 0:
V_h = \left\{ v_h \in X_h : v_h(\Gamma_D) = 0 \right \}
The basis of V_h can be expressed as
\varphi_h \in X_h : \varphi_j(x_i) = \delta_{ij} = \begin{cases} 0 & i \ne j \\ 1 & \text{otherwise} \end{cases}
Therefore
v_h(x) = \sum_{j=1}^{N_h} v_j \varphi_j(x) \implies v_h(x_i) = v_i
The error in the discretization is bounded: \|u - u_h\| \le C \cdot h.
To be continued
Multidimensional elliptic PDEs
Multidimensional PDEs are just like monodimensional PDEs, but multidimensional:
\begin{cases} Lu := -\operatorname{div}(\mu \nabla u) + \vec{b} \nabla u + \sigma v \\ u = 0 & \text{on } \Gamma_D \\ \mu \frac{\partial u}{\partial n} = \phi & \text{on } \Gamma_N \end{cases}
Integrals become multidimensional, derivatives become gradients and the integration by parts is performed through the Green formula:
-\int_{\Omega} \operatorname{div}(\mu \nabla u) v= \int_{\Omega} mu \nabla u \nabla v - \int_{\partial\Omega} \mu \underbrace{\nabla u \vec{n}}_{\frac{\partial u}{\partial n}} v
Massaging the problem (if you get this joke, you are authorised to open the source of this file and to add your name in the comment below this paragraph) as in the monodimensional case and using the Green formula instead of the integration by parts, we can write that
\begin{align*} a(u, v) &= \int_{\Omega} \mu \nabla u \nabla v \int_{\Omega} \vec{b} \nabla u v + \int_{\Omega} \sigma u v \\ F(v) &= \int_{\Omega} fv + \int_{\partial\Omega} \mu \frac{\partial u}{\partial n} v \\ &= \int_{\Omega} fv + \underbrace{\int_{\Gamma_D} \mu \frac{\partial u}{\partial n} v}_{=0} + \int_{\Gamma_N} \underbrace{\mu \frac{\partial u}{\partial n}}_{\psi} v \\ &= \int_{\Omega} fv + \int_{\Gamma_N} \psi v \end{align*}
Continuing with the massage we can get to the weak form, to the Galerkin approximation and finally, to the algebraic system.
The construction of V_h follows the same logic: the definitions are the same as in the monodimensional case with \Omega as the domain. It is notable that X_h \sub H'(\Omega) and V_h \sub H'_{\Gamma_D} where
V \equiv H^1_{\Gamma_D}(\Omega) = \left\{ v \in H'(\Omega), v|_{\Gamma_D} = 0 \right\} \\ H^1(\Omega) = \left\{ v : \Omega \to \mathbb{R}, v \in L^2(\Omega), \frac{\partial v}{\partial x_j} \in L^2(\Omega) \ \forall j = 1, 2, \dots, n \right\}
Lax-Milgram lemma
Assume that V is an Hilbert space with norm \|\cdot\|, F is a linear functional on V (i.e. F \in V' and bounded), a is a bilinear, continuous and coercive form (i.e. \exists u \gt 0 : |a(u, v)| \le u \|u\|\|v\|\ \forall u, v \in V, \exists \alpha \gt 0 : a(v, v) \ge \alpha \|v\|^2\ \forall v \in V).
Then there exists a unique solution u to the weak problem that is also bounded (i.e. \|u\| \lt \frac{1}{\alpha}\|F\|_{V'})
In the previous lemma, V' is the dual space of V and is defined as
V' = \left\{ F : V \to \mathbb{R} : F \text{ is linear, bounded and } \|F\|_{V'} = \frac{|F(v)|}{\|v\|_V} \le C \right\}
If the assumptions of the Lax-Milgram lemma are satosfied, then the Galerkin problem has a unique solution and is bounded independently of h (the solution is stable):
\|u_h\| \lt \frac{1}{\alpha}\|F\|_{V'}
Moreover, the resulting linear alegraic system will be nonsingular.
We will not proove the Lax-Milgram lemma but we will only see the conditions a specific problem must satisfy to satisfy the assumptions of the lemma.
Assume we have to proove that the Lax-Milgram lemma applies to this problem (homogeneous mixed problem):
a(u, v) = \int_a^b \mu u' v' + b u' v + \sigma u v \\ F(v) = \int fv + \gamma v(1) \\ \Omega = {x \in \mathbb{R} : a \le x \le b}, a = 0, b = 1
The Lax-Milgram lemma requires that
- V is an Hilbert space
- a is bilinear
- a is continuous
- a is coercive
- F is linear
- F is bounded
- f \in L^2(\Omega)
For what concerns point (1), we take for granted that V is an Hilbert space. We have to choose one of two norms to use in the proofs. We can either choose the complete norm
\|v\| = \sqrt{\|v\|^2_{L^2}(0, 1) + \|v'\|_{L^2(0, 1)}}
or the reduced norm
|v| = \sqrt{\|v'\|_{L^2(0, 1)}}
The reduced norm can be used only with a non-empty Dirichlet boundary. For this proof, we chose to use the complete norm.
Proving (2), (5) and (7) is trivial.
We will now prove (3). By definition, a is continuous if \exists M \gt 0 : |a(u, v)| \le M \|u\| \|v\| \ \forall u, v \in V.
|a(u, v)| \le \left| \int_a^b \mu u' v' \right| + \left| \int_a^b b u' v \right| + \left| \int_a^b \sigma u v \right| = (*)
We can write that
\left| \int_a^b \mu u' v' \right| \le \underbrace{\max_{x} |\mu(x)|}_{\mu_{max}} \int_a^b |u' v'| \le \mu_{max} \|u'\|_{L^2(a, b)} \|v'\|_{L^2(a, b)} \\ \left| \int_a^b b u' v \right| \le \underbrace{\max_x |b(x)|}_{b_{max}} \int_a^b |u'v| \le b_{max} \|u'\|_{L^2(a, b)} \|v'\|_{L^2(a, b)} \\ \left| \int_a^b \sigma u v \right| \le \underbrace{\max_x |\sigma(x)|}_{\sigma_{max}} \int_a^b |u v| \le \sigma_{max} \|u\|_{L^2(a, b)} \|v\|_{L^2(a, b)}
Knowing that \|v\|_{L^2} \lt \|v\| and that \|v'\|_{L^2} \lt \|v\|, we can conclude that
\begin{align*} (*) &\le \mu_{max} \|u'\|_{L^2(a, b)} \|v'\|_{L^2(a, b)} + b_{max} \|u'\|_{L^2(a, b)} \|v'\|_{L^2(a, b)} + \sigma_{max} \|u\|_{L^2(a, b)} \|v\|_{L^2(a, b)} \\ &\le \mu_{max} \|u\| \|v\| + b_{max} \|u\| \|v\| + \sigma_{max} \|u\| \|v\| \\ &= \underbrace{(\mu_{max} + b_{max} + \sigma_{max})}_{M} \|u\| \|v\| \\ &= M \|u\| \|v\| \end{align*}
Assuming that \mu, b and \sigma are bounded from above, we can say that a is continuous.
We will now prove (4). By definition, a is coercive if \exists \alpha : a(v, v) \ge \alpha \|v\|^2.
We know that a(v, v) = \int_a^b \mu(v')^2 + \int_a^b b v' v + \int_a^b \sigma v^2.
Assume that \mu \gt 0 then
\int_a^b \mu (v')^2 \ge \mu_{min} \int_a^b (v')^2 = \mu_{min} \|v'\|_{L^2(a, b)}^2
Assume that b(1) \gt 0 then
\int_a^b b v' v = \int_a^b b \frac{1}{2} (v^2)' = \frac{1}{2} \int_a^b b (v^2)' = -\frac{1}{2} \int_a^b b v^2 + \underbrace{\frac{1}{2} [b v^2]_a^b}_{\gt 0} \ge -\frac{1}{2} \int_a^b b' v^2
Assume \delta = \sigma - \frac{1}{2}b' \ge 0 then
\begin{align*} a(v, v) &= \mu_{min} \|v'\|_{L^2(a, b)}^2 \int_a^b (\sigma - \frac{1}{2}b') v^2 \\ &\ge \mu_{min} \|v'\|_{L^2(a, b)} + \delta \int_a^b v^2 \\ &= \mu_{min} \|v'\|_{L^2(a, b)}^2 + \delta \|v\|_{L^2(a, b)}^2 \\ &\ge \underbrace{\min(\mu_{min}, \delta)}_{\alpha \gt 0} \cdot (\|v'\|_{L^2(a, b)} + \|v\|_{L^2(a, b)}^2) \\ &= \alpha\|v\|^2 \end{align*}
We will now prove (6). Since \|1\|_{L^2(a, b)} = \sqrt{\int_a^b 1^2 dx} = \sqrt{b - a} then
v(b) = v(a) + \int_a^b v'(x) dx \implies |v(b)| = \cancel{|v(a)|} \int_a^b v'(x) \cdot 1 dx \le \|v'\|_{L^2(a, b)} \|1\|_{L^2(a, b)} = \|v'\|_{L^2(a, b)} \sqrt{b - a} \le \sqrt{b - a} \|v\|
this means that
\begin{align*} |F(v)| &= \left|\int_a^b fv + \gamma v(b) \right| \\ &\le \left| \int_a^b fv \right| + |\gamma v(1)| \\ &\le \|f\|_{L^2(a, b)} \|v\|_{L^2(a, b)} + |\gamma| |v(b)| \\ &\le (\|f\|_{L^2(a, b)} + |\gamma|) \|v\| \\ &= \alpha \|v\| \end{align*}
We can now conclude that
\|F\|_{V'} = \sup_{v \in V} \frac{F(v)}{\|v\|} \le \alpha
Since f and \gamma are given, the problem is said to be bounded by the data.
Since the assumptions of the Lax-Milgram lemma are satisfied, then we can say that there exists a unique solution for a(u, v) = F(v) \ \forall v \in V that is also bounded.
The Lax-Milgram lemma also holds with multidimensional domains.
If we take u = v then
\alpha \|u\|^2 \le a(u, u) = F(u) = \int fu + \gamma u(b) \le \|f\|_{L^2(a, b)} \|u\|_{L^2(a, b)} + |\gamma| \|u\| \le (\|f\|_{L^2(a, b)} + |\gamma|)\|u\|
therefore
\|u\| \le \frac{1}{\alpha} (\|F\|_{L^2(a, b)} + |\gamma|)
This means that the solution is bounded independently on the value of h (stability property).
TODO: NDim generalization
Non homogeneous elliptic PDEs
Up until now, we have assumed that the Dirichlet boundary conditions were shaped like u(\Gamma_D) = 0. This is called homogeneous problem. If u(\Gamma_D) = \varphi then it is called heterogeneous problem.
In the case of an heterogeneous problem, we introduce a lifting R of \varphi such that u^0 = u - R\varphi and u^0(a) = 0.
The problem now become
\begin{align*} a(u, v) &= F(v) \\ a(u^0 + R\varphi, v) &= F(v) \\ a(u^0, v) + \underbrace{a(R\varphi, v)}_{G(v)} = F(v) \\ a(u^0, v) = F(v) - G(v) \end{align*}
A good choice for R in 1D is a choice s.t. R\varphi = \varphi. In the multidimensional case, we choose an R such that u(\Gamma_D) = 0.
Galerkin approximation
TODO
Practical numerical solution
TODO
Parabolic PDEs
Parabolic PDEs are a generalization of elliptic PDEs that differs only for the fact that they are time-dependent.
A general parabolic problem can be expressed as
\text{$\forall t \in (0, T)$, find $u$ : }\begin{cases} \frac{\partial u}{\partial t} + Lu = f & \text{in } \Omega \\ \text{Boundary conditions} & \text{in } \partial\Omega, 0 \lt t \lt T \\ \text{Initial condition} & \text{in } x \in \Omega, t = 0 \end{cases}
where L is the elliptic operator already defined for elliptic PDEs and so are the boundary conditions, except for the fact that they must be true \forall t. The solution u also depends on time.
The time-dependent domain where the PArabolid PDE applies is defined as Q = \Omega \times (0, T) and, technically, it is a cylinder.
The procedure to get to “something that can be computed” is really similar to the one used for elliptic PDEs.
flowchart LR
id1[Strong parabolic problem]
id2[Weak parabolic problem]
id3[Galerkin problem]
id4[Algebraic problem]
id1-->id2-->id3-->id4In order to get to the weak parabolic problem we first multiply by a test function v and then integrate, moving all the known terms to the right hand side of the equation.
\frac{\partial u}{\partial t} -\operatorname{div}(\mu \nabla u) + \vec{b} \nabla u + \sigma u = f \\ \int_\Omega \frac{\partial u}{\partial t} v \ dx - \int_\Omega \operatorname{div}(\mu \nabla u) v \ dx + \int_\Omega \vec{b} \nabla u v \ dx + \int_\Omega \sigma uv = \int_\Omega fv \ dx \qquad \forall v \in V\\ \int_\Omega \frac{\partial u}{\partial t} v \ dx + \underbrace{\int_\Omega \mu \nabla u \nabla v + \int_\Omega \vec{b} \nabla u v \ dx + \int_\Omega \sigma u v \ dx}_{a(u, v)} = \underbrace{\int_\Omega f v \ dx + \int_{\Gamma_N} \psi v \ d\gamma}_{F(v)} \qquad \forall v \in V \\ \int_\Omega \frac{\partial u}{\partial t} v \ dx + a(u, v) = F(v) \qquad \forall v \in V
Therefore, the weak parabolic form of the problem is
\text{$\forall t \in (0, T)$, find $u \in V$ s.t. } \int_\Omega \frac{\partial u}{\partial t} v \ dx + a(u, v) = F(v) \qquad \forall v \in V
where V is defined exactly in the same way it was defined for elliptic PDEs.
Let V_h \sub V, \dim(V_h) = N_h then we can get to the Galerkin approximation of the problem:
\text{$\forall t \in (0, T)$, find $u_h \in V_h$ s.t. } \begin{cases} \int_\Omega \frac{\partial u_h}{\partial t} v_h \ dx + a(u_h, v_h) = F(v_h) \quad \forall v_h \in V_h \\ u_n|_{t=0} = u_{0h} \in V_h \end{cases}
Knowing that
v_h(x) = \sum_{j = 1}^{N_h} v_j \varphi_j(x)
we can apply even more algebraic transformations:
\int_\Omega \frac{\partial}{\partial t} \left( \sum_{j=1}^{N_h} u_j(t) \varphi_j \right) \varphi_i \ dx + a\left( \sum_{j = 1}^{N_h} u_j(t \varphi_j, \varphi_i) \right) = F(\varphi_i) \qquad \forall i = 1, 2, \dots, N_h \\ \sum_{j=1}^{N_h} \frac{\partial}{\partial t} u_j(t) \int_\Omega \varphi_j \varphi_i \ dx + \sum_{j=1}^{N_h} u_j(t) a(\varphi_j, \varphi_i) = F(\varphi_i) \qquad \forall i = 1, 2, \dots, N_h \\ M \frac{\partial \vec{u}}{\partial t} + A \vec{u} = \vec{F}
In the last step we got to a system of ordinary differential equations (the actual Galerkin approximation).
The Galerkin approximation is also called semidiscretization of the weak parabolic problem because the time is still continuous.
To get to the fully discretized problem, the \theta-method must be applied.
Let y^n = y(t^n) then \frac{\partial y}{\partial t} in an ODE can be approximated as follows:
\frac{y^{n+1} - y^{n}}{\Delta t} = \theta f(t^{n+1}, y^{n+1}) + (1 - \theta)f(t^n, y^n) \qquad 0 \le \theta \le 1
If \theta = 0 we get the forward Euler method, if \theta = \frac{1}{2} we get the Crank-Nicolson method and if \theta = 1 we get the backward Euler method. The Crank-Nicolson method gives a second order convergence, while the Euler methods only gives first order convergence.
With the application of the \theta-method, we can finally write the fully discretized problem:
M \frac{\vec{u}^{n+1} - \vec{u}^n}{\Delta t} + \theta A \vec{u}^{n+1} + (1 - \theta) A \vec{u}^n = \theta \vec{F}^{n+1} + (1 - \theta)\vec{F}^{n} \qquad n = 1, 2, \dots, N \\ \underbrace{\left( \frac{1}{\Delta t}M + \theta A \right)}_{K} \vec{u}^{n+1} = \underbrace{\theta \vec{F}^{n+1} + (1 - \theta) \vec{F}^n + \frac{1}{\Delta t} M \vec{u}^n - (1 - \theta) A \vec{u}^n}_{\vec{G}^n} \qquad n = 1, 2, \dots, N \\ K \vec{u}^{n+1} = \vec{G}^n
The fully discretized form is composed by a linear system for each timestep n. Since the solution of the n+1th timestep depends on the solution of the nth one, the linear systems must be resolved sequentially.
As K does not depend on time, it may be useful to apply LU or Cholesky factorization on it once and for all before solving the system for each timestep. For a big enough number of timesteps, the time spent performing the factorization will be masked by the performance gain.
Well-posedness of parabolic PDE problems
For the well-posedness of elliptic problems, coercivity was required. For parabolic problems, weak coercivity is enough.
Assume we want to prove the well-posedness of the following problem
\begin{cases} \frac{\partial u}{\partial t} - \nu \frac{\partial^2 u}{\partial x^2} = f & \bar a \lt x \lt \bar b \\ \nu \frac{\partial u}{\partial x}(\bar a) = \frac{\partial u}{\partial x}(\bar b) = 0 & \Gamma_N \\ u(x, 0) = u_0(x) & \bar a \lt x \lt \bar b \end{cases}
with \nu \gt 0 and constant.
If we multiply by a test function and integrate, we get that
\int_{\bar a}^{\bar b} \frac{\partial u}{\partial t} v + \underbrace{\int_{\bar a}^{\bar b} \nu \frac{\partial u}{\partial x} \frac{\partial v}{\partial x}}_{a(u, v)} - \underbrace{\left[ \nu \frac{\partial u}{\partial x} v \right]_{\bar a}^{\bar b}}_{0} = \int_{\bar a}^{\bar b} fv
Consider a(v, v) = \int_{\bar a}^{\bar b} \nu (\frac{\partial v}{\partial x})^2 = \nu \|\frac{\partial v}{\partial x}\|_{L^2(\bar a, \bar b)}^2: there cannot be coercivity because a(v, v) cannot be, in any way, larger than the complete norm of v squared.
Remember that the complete norm of v is expressed as
\|v\| = \sqrt{\|v\|_{L^2(\bar a, \bar b)}^2 + \left\| \frac{\partial v}{\partial x} \right\|_{L^2(\bar a, \bar b)}^2}
Consider instead a(v, v) + \lambda \|v\|_{L^2(\bar a, \bar b)}^2 then
a(v, v) + \lambda \|v\|_{L^2(\bar a, \bar b)} = \nu \left\| \frac{\partial u}{\partial x} \right\|_{L^2(\bar a, \bar b)}^2 + \lambda \|v\|_{L^2(\bar a, \bar b)}^2 \ge \underbrace{\min(\nu, \lambda)}_\alpha \cdot \underbrace{\left( \|v\|_{L^2(\bar a, \bar b)}^2 + \left\| \frac{\partial v}{\partial x}_{L^2(\bar a, \bar b)}^2 \right\| \right)}_{\|v\|^2} = \alpha \|v\|^2
hence, we got that a is weakly coercive.
Stability of the \theta-method
The term stability refers to the fact that the solution does not explode and that is bounded in time.
Depending on the value of \theta, we can either have unconditional stability (the method is stable for an arbitrary \Delta t) or conditional stability (the method is stable only for small enough values of \Delta t).
To study the stability of the fully discretized problem, we apply the \theta-method to the Galerkin problem:
\left( \frac{u_h^{k+1} - u_h^{k}}{\Delta t}, v_h \right) + a(\theta u_h^{k+1} + (1 - \theta)u_h^k, v_h) = \theta F^{k+1}(v_h) + (1 - \theta)F^k(v_h)
for each timestep k and where (\cdot, \cdot) is the sclar product:
(a, b) = \int a(x) \cdot b(x) \ dx
Assume for simplicity that F = 0 and that \theta = 1 (backwrd Euler). As we want to use coercivity, we take a(u_h^{k+1}, u_h^{k+1}).
(u_h^{k+1}, u_h^{k+1}) + \Delta t \underbrace{a(u_h^{k+1}, u_h^{k+1})}_{\ge \alpha \|u_h^{k+1}\|_V^2} = (u_h^k, u_h^{k+1}) \\ \|u_h^{k+1}\|_{L^2(\Omega)}^2 + \Delta t \alpha \|u_h^{k+1}\|_V^2 \le \|u_h^k\|_{L^2(\Omega)} \|u_h^{k+1}\|_{L^2(\Omega)} \le \frac{1}{2} \|u_h^k\|^2_{L^2(\Omega)} + \frac{1}{2} \|u_h^{k+1}\|_{L^2(\Omega)}^2
where we used the coercivity of a (even if it is weakly coercive, it works anyway).
We will now prove boundedess:
\|u_h^{k+1}\|_{L^2(\Omega)}^2 + 2 \Delta t \alpha \|u_h^{k+1}\|_{L^2(\Omega)}^2 \le \|u_h^k\|_{L^2(\Omega)}^2 \\ \left( \|u_h^{k+1}\|_{L^2(\Omega)}^2 - \|u_h^k\|_{L^2(\Omega)}^2 \right) + 2 \Delta t \alpha \|u_h^{k+1}\|_{L^2(\Omega)}^2 \le 0 \\ \left( \|u_h^{n+1}\|_{L^2(\Omega)}^2 - \|u_h^0\|_{L^2(\Omega)}^2 \right) + 2 \alpha \Delta t \sum_{k = 0}^n \|u_h^{k+1}\|_V^2 \le 0 \\ \underbrace{\|u_h^{n+1}\|_{L^2(\Omega)}^2 + 2 \alpha \Delta t \sum \|u_n^{k+1}\|_{L^2(\Omega)}^2}_{\text{One possible norm for $u_h^{n+1}$}} \le \|u^0_h\|_{L^2(\Omega)}
where we used the fact that the V-norm of somethins is smaller or equal to the L^2 norm of the same thing.
We will now proove stability
(1 + 2 \alpha \Delta t) \|u_h^{k+1}\|_{L^2(\Omega)}^2 \le \|u_h^k\| \\ \|u_h^k\|_{L^2(\Omega)} \le \left( \frac{1}{\sqrt{1 + 2 \alpha \Delta t}} \right)^k \|u_h^0\|_{L^2(\Omega)}
therefore the backward Euler method is unconditionally stable.
In the previous example we used the fact that a even if it was declared as weakly coercive. This is perfectly fine because if we perform a specific change of variables on a weakly coercive form, we get a normally coercive form:
Take a weakly coercive form a:
\frac{\partial u}{\partial t} - \nu \frac{\partial^2 u}{\partial x^2} = f \implies a(u, v) = \int_\Omega \frac{\partial u}{\partial x} \frac{\partial v}{\partial x}
If we take u = e^{\lambda t}w then the problem become
\frac{\partial}{\partial t} \left( e^{\lambda t} w \right) - \nu \frac{\partial^2}{\partial t^2} \left( e^{\lambda t} w \right) = f \implies \tilde a(w, v) = \int_{\Omega} \frac{\partial w}{\partial x} \frac{\partial v}{\partial x} + \lambda \int_{\Omega} w^2 = a(w, v) + \lambda \|v\|_{L^2(\Omega)}^2 \ge \lambda \|w\|_{L^2(\Omega)}^2
While the backward Euler is unconditionally stable, in general, this property depends on the chosen value for \theta, in particular, the \theta-method is unconditionally stable only for \theta \ge 0, otherwise we must have that \Delta t is small enough:
\Delta t \lt \frac{2}{(1 - 2\theta) \lambda_{max}(A)} \simeq \frac{2}{1 - 2\theta} h^2
Let A and M be two matrices. A generalized eigenvalue problem consists in finding all the eigenpairs (\lambda, v) s.t. Av = \lambda M v. If A is SPD then all the eigenvalues \lambda are strictly positive and real numbers.
Let a \in V \times V be a bilinear form. The eigenpairs of a can be written as all the tuples (\lambda, w) s.t.
a(w, v) = \lambda(w, v) = \lambda \int_{\Omega} w v \qquad \forall v \in V
If, instead of V, we take V_h, we get that (\lambda_h, w_h) is an approzimation of (\lambda, w) if
a(w_h, v_h) = \lambda_h(w_h, v_h) \qquad \forall v_h \in V_h \\ a\left( \sum_j w_j \varphi_j, \varphi_i \right) = \lambda \left( \sum_j w_j \varphi_j, \varphi_i \right) \qquad \forall i = 1, 2, \dots, N_h\\ \sum_j w_j a(\varphi_j, \varphi_i) = \lambda_h \sum_j w_j (\varphi_j, \varphi_i) \qquad \forall i = 1, 2, \dots, N_h \\ A\vec{w} = \lambda_h M \vec{w}
therefore the problem of the approximation of the eigenpairs of a bilinear form is equivalent to a generalized eigenvalue problem.
Let (w_h^j, \lambda_h^j) be the jth eigenpair and assume that w_h^k is normalized. It holds that
(w_h^j, w_h^i) = \delta_{ij} \qquad \|w_h^j\|_{L^2(\Omega)}^2 = 1
Let u_h be the approximation of the solution, v_h be the test function and F = 0, then, if we apply the \theta-method, we get that
\int_\Omega \frac{u_h^{k+1} - u_h^k}{\Delta t} v_h + a(\theta u_h^{k+1} + (1 - \theta) u_h^k, v_h) = 0 \qquad \forall v_h \in V_h \\ \int_\Omega \frac{1}{\Delta t} \left( \sum_j (u_j^{k+1} - u_j^k) w_h^j \right) w_h^i + a\left(\sum_j u_j^{k+1} w_h^j + (1 - \theta)\sum_j u_j^k w_h^j, w_h^i\right) = 0 \qquad \forall i = 1, 2, \dots, N_h \\ \frac{1}{\Delta t} \sum_j (u_j^{k+1} - u_j^k) \underbrace{\int_\Omega w_h^j, w_h^i}_{\delta_{ij}} + \sum_j (\theta u_j^{k+1} - (1 - \theta) u_j^k) \underbrace{a(w_h^j, w_h^i)}_{\lambda_h^j(w_h^j, w_h^i) = \lambda_h^j \delta_{ij}} = 0 \qquad \forall i = 1, 2, \dots, N_h \\ \frac{1}{\Delta t} (u_i^{k+1} - u_i^{k}) + (\theta u_i^{k+1} + (1 - \theta) u_i^k) \lambda_h^i = 0 \qquad \forall i = 1, 2, \dots, N_h \\ u_i^{k+1} + \Delta t \theta u_i^{k+1} \lambda_h^i = u_i^k - \Delta t (1 - \theta) \lambda_h^i u_i^k \qquad \forall i = 1, 2, \dots, N_h \\ (1 + \Delta t \theta \lambda_k^i) u_i^{k+1} = [1 - \Delta t (1 - \theta) \lambda_h^i] u_i^k \qquad \forall i = 1, 2, \dots, N_h \\ u_i^{k+1} = \underbrace{\frac{1 - \Delta t (1 - \theta) \lambda_h^i}{1 + \Delta t \theta \lambda_h^i}}_{\rho_i} u_h^k \qquad \forall i = 1, 2, \dots, N_h
Assume |\rho_i| \lt 1 then |u_i^{k+1}| \lt |u_i^k| \lt \dots \lt |u_i^0|, therefore
\sum_{i = 1}^{N_h} |u_i^n|^2 \lt \sum_{i = 1}^{N_h} |u_i^0|^2
Since
u_h^n(x) = \sum_{j = 1}^{N_h} u_j^n w_h^j(x)
then
\|u_h^n\|_{L^2(\Omega)}^2 = \int_\Omega (u_h^n)^2 = \int_\Omega \left( \sum u_j^n w_h^j \right)\left( \sum u_j^n w_h^j \right) = \sum_j \sum_i u_j^n u_i^n \underbrace{\int_\Omega w_h^j w_h^i}_{\delta_{ij}} = \sum_j (u_j^n)^2
therefore
\|u_h^n\|_{L^2(\Omega)} \le \|u_h^n\|_{L^2(\Omega)}
hence, stability of the \theta-method.
To prove stability, we assumed that |\rho_i| < 1, we will now study when that statement is true.
|\rho_i| \lt 1 \iff -1 - \theta \lambda_h^i \Delta t \lt 1 - (1 - \theta) \lambda_h^i \Delta t \lt 1 + \theta \lambda_h^i \Delta t \iff 2 \theta - 1 \gt - \frac{2}{\lambda_h^i \Delta t}
If \theta \ge \frac{1}{2} the inequality is always satisfied and the method in unconditionally stable, otherwise we must impose that
\Delta t \lt \frac{2}{(1 - 2 \theta) \lambda_h^i}
Since that condition should hold for all the eigenvalies \lambda_h^i, it will suffice to impose it only for the bigger one, so
\Delta t \gt \frac{2}{(1 - 2 \theta) \lambda_h^{N_h}} \simeq \frac{2}{1 - 2\theta} h^2
Domain decomposition
Domain Decomposition methods are used to split a big problem into two or more smaller problems. Assume you have a problem so solve that looks like
\begin{cases} Lu = f & \Omega = \{x : \bar a \lt x \lt \bar b\}\\ u(\bar a) = \gamma_1 \\ u(\bar b) = \gamma_2 \end{cases}
We want to split \Omega into two subdomains so that we can get a new, equivalent problem:
\begin{cases} Lu_1 = f & \Omega_1 = \{x : \bar a \lt x \lt \Gamma\} \\ Lu_2 = f & \Omega_2 = \{x : \Gamma \lt x \lt \bar b\} \\ u_1(\bar a) = \gamma_1 \\ u_2(\bar b) = \gamma_2 \\ u_1(\Gamma) = u_2(\Gamma) \\ (\mu u_1')(\Gamma) = (\mu u_2')(\Gamma) \end{cases}
\Gamma is called the interface between \Omega_1 and \Omega_2.
We see in the new, equivalent problem, that on the interface, we need to impose continuity of both the solution and its derivative (interface conditions). Interface conditions are both Dirichlet and Neumann.
If we define
u = \begin{cases} u_1 & \Omega_1 \\ u_2 & \Omega_2 \end{cases}
then we can recover the solution of the original problem from the unknowns of the equivalent one.
u_1 and u_2 are called local solutions.
We use an iterative method (the Dirichlet-Neumann DD method) to solve the new equivalent problem splitting it into two separate problems:
\begin{cases} Lu_1^{(k)} = f & \Omega_1 \\ u_1^{(k)} = u_2^{(k-1)} & \Gamma \\ \text{Boundary conditions} & \partial\Gamma_1 \cap \partial\Gamma \end{cases} \\ \begin{cases} Lu_2^{(k)} = f & \Omega_2 \\ \mu u_2'^{(k)} = \mu u_1'^{(k)} & \Gamma \\ \text{Boundary conditions} & \partial\Gamma_2 \cap \partial\Gamma \end{cases}
Lets call P_1 the problem on \Omega_1 and P_2 the one on \Omega_2.
The DN-DD method is sequential: u_2^{(k)} (from P_2) depends on u_1^{(k)} (from P_1) which in turn depends on u_2^{(k-1)} (from P_2 on the previous iteration) and so on.
We modify P_2 to make the two P_1 and P_2 independent so that they can be solved in parallel (in red the changes to the algorithm):
\begin{cases} Lu_2^{(k)} = f & \Omega_2 \\ \mu u_2'^{(k)} = \mu \textcolor{red}{u_1'^{(k - 1)}} & \Gamma \\ \text{Boundary conditions} & \partial\Gamma_2 \cap \partial\Gamma \end{cases}
With k \to \infty, we have that u_i^{(k)} \to u_i.
Depending on the choice for \Gamma, the algorithm may or may not converge. To improve convergence, we modify P_1 to have a relaxation parameter (changes in red) to get the DN-\theta problem:
\begin{cases} Lu_1^{(k)} = f & \Omega_1 \\ u_1^{(k)} = \textcolor{red}{\theta u_2^{(k-1)} + (1-\theta) u_1^{(k-1)}} & \Gamma \\ \text{Boundary conditions} & \partial\Gamma_1 \cap \partial\Gamma \end{cases}
where \theta > 0 is called relaxation parameter and can be chosen arbitrarily.
With 0 \lt \theta \lt \theta_{max}, convergence is guaranteed.
With a 1D problem, splitting it in two domains, it exists \theta_{opt} s.t. the algorithm converges in one single iteration.
Appendix
Normal derivative
Let v : \Omega \sub \mathbb{R}^d \to \mathbb{R}, then the normal derivative of v is defined as
\frac{\partial v}{\partial n} = \nabla v \cdot n
where n is the normal direction of v in each point.
Divergence
Let \vec{w} \in \mathbb{R}^d, d \in \mathbb{N}^+ then the divergence operator applied to \vec{w} is defined as
\operatorname{div}(\vec{w}) = \sum_{i = 1}^d \frac{\partial w_i}{\partial x_i}
Gradient
Let v : \Omega \sub \mathbb{R}^d \to \mathbb{R}, d \in\mathbb{N}^+, then the gradient of v is defined as
\nabla v = \begin{bmatrix} \frac{\partial v}{\partial x_1} \\ \frac{\partial v}{\partial x_2} \\ \vdots \\ \frac{\partial v}{\partial x_d} \\ \end{bmatrix} \in \mathbb{R}^d
Let \vec{x} \in \Omega then \nabla v(\vec{x}) gives the direction of steepest ascent. If \nabla v(\vec{x}) = 0 then \vec{x} can be either a local maximum, a local minimum or a saddle point.
Laplacian
Let v : \Omega \sub \mathbb{R}^d \to \mathbb{R}, d \in \mathbb{N}^+, then the laplacian of v is defined as
\Delta v = \sum_{i = 1}^d \frac{\partial^2 v}{\partial x_i^2}
Distributional derivative
Let v : \mathbb{R}^n \to \mathbb{R}, then D is the distributional derivative of v if
\int_{-\infty}^{+\infty} Dv \cdot w \, dx = -\int_{-\infty}^{+\infty} v \frac{dw}{dx} \, dx \qquad \forall w \in \mathcal{C}^\infty(\mathbb{R}) \text{ s.t. } w : \mathbb{R} \to \mathbb{R}, \lim_{x \to \pm \infty}w(x) = 0
If v \in \mathcal{C}^1 then the distributional derivative is the same as the conventional derivative.
Intuitively, in the case of jump discontinuities, the distributional derivative is a piecewise derivative where at each discontinuity point \bar x, we can choose either v(\bar x) = v(\bar x^-) or v(\bar x) = v(\bar x^+).
Distributional derivatives are also called weak derivatives.
Square-integrable functions
The space of square-integrable functions L^2 on \Omega \sub \mathbb{R}^n is defined as
L^2(\Omega) = \left\{ f : \Omega \to \mathbb{R} \middle| \int_\Omega f(x)^2 \, d\Omega \lt +\infty \right\}
The integral is to be considered a Lebesgue integral.
If u, v \in L^2 then u' \cdot v' \in L^1.
L^2 is an Hilbert space where
\lang f, g \rang_{L^2(\Omega)} = \int_\Omega f(x)g(x) \, d\Omega \\ \|f\|_{L^2(\Omega)} = \sqrt{\lang f, f \rang_{L^2(\Omega)}}
Forms
A form a is a relation a : V \times V \to \mathbb R (where V is a function space).
A form a is called bilinear if
a(\lambda u + \mu w, v) = \lambda a(u, v) + \mu a(w, v) \\ a(u, \lambda w + \mu v) = \lambda a(u, w) + \mu a(u, v)
A form a is called coercive if
\exists \alpha \gt 0 : a(v, v) \gt \alpha \|v\|_V^2
A form a is called weakly coercive if
\exists \lambda \ge 0, \alpha \gt 0 : a(v, v) + \lambda \|v\|_{L^2}^2 \ge \alpha \|v\|_V^2
Functionals
TODO
Cauchy-Schwarz inequality
Let u, v \in V, then
\left|\int_\Omega uv\right| \le \|u\|_{L^2(\Omega)} \|v\|_{L^2(\Omega)} \\ \left|\int_\Omega u'v'\right| \le \|u'\|_{L^2(\Omega)} \|v'\|_{L^2(\Omega)} \\ \left|\int_\Omega u'v\right| \le \|u'\|_{L^2(\Omega)} \|v'\|_{L^2(\Omega)} \\
Young inequality
Let A, B \in \mathbb{R} then
AB \le \varepsilon A^2 + \frac{1}{4\varepsilon}B^2 \qquad \forall \varepsilon \gt 0
In particular, for \varepsilon = \frac{1}{2}, then
AB \le \frac{A^2}{2} + \frac{B^2}{2}